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ABSTRACT 

Wc investigate the linear stability of a shocked accretion flow onto a black hole 
in the adiabatic limit. Our linear analyses and numerical calculations show 
that, despite the post-shock deceleration, the shock is generally unstable to 
non-axisymmetric perturbations. The simulation results of Moltcni, Toth & 
Kuznetsov can be well explained by our linear eigenmodes. The mechanism 
of this instability is confirmed to be based on the cycle of acoustic waves 
between the corotation radius and the shock. We obtain an analytical formula 
to calculate the oscillation period from the physical parameters of the flow. 
We argue that the quasi-periodic oscillation should be a common phenomenon 
in accretion flows with angular momentum. 

Key words: accretion, accretion discs - black hole physics - hydrodynamics 
- instabilities - shock waves 



1 INTRODUCTION 

Hydrodynamic instabilities of shocked accretion flows may explain quasi-periodic oscilla- 
tion (QPO) processes occurring in black hole candidates. The structure of stationary black 
hole accretion flows involving standing shocks was first described by Fukue (1987). Subse- 
quently, shock studies have been made extensively in both inviscid and viscous accretion 
flows (Chakrabarti & Das 2004; Gu & Lu 2004 and references therein). However, even with 
the simple inviscid hypothesis, the stability of the shock is not fully understood. In the 
isothermal limit, Nakayama (1992) introduced a global instability between a sonic point and 
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a shock, and found the criterion that "post-shock acceleration causes instability", which 
was confirmed by the simulations of Nobuta & Hanawa (1994). Moreover, Nakayama (1994) 
investigated this instability in an adiabatic flow and claimed that such a criterion is also 
correct unless the shock is extremely strong. All the above works, however, were only for 
axisymmetric perturbations. 

The pioneering work of Papaloizou & Pringle (1984) found a non-axisymmetric insta- 
bility based on the acoustic cycle between the corotation radius and the boundary. The 
Papaloizou-Pringle instability (hereafter PPI) is known to take place in discs or tori, in 
which the radial velocity is initially zero. The mechanism of such an instability has been dis- 
cussed extensively by Goldreich & Narayan (1985), Goldreich, Goodman & Narayan (1986), 
Narayan, Goldreich & Goodman (1987) and Kato (1987). The effect of radial advection on 
the PPI was investigated by Blaes (1987), who found that the PPI is strongly stabilized by 
advection at the inner boundary. Another type of non-axisymmetric instability was found 
in a spherical accretion flow by Foglizzo & Tagger (2000) and Foglizzo (2001, 2002), which 
is based on the cycle of entropy/vorticity perturbations and acoustic waves in the subsonic 
region between a stationary shock and a sonic surface. Recently, Gu & Foglizzo (2003, here- 
after Paper I) studied a shocked accretion disc in the isothermal limit and found that the 
shock is generally linearly unstable to non-axisymmetric perturbations despite the post- 
shock deceleration. Paper I pointed out that such an instability is a form of PPI modified by 
advection and the presence of the shock. Apart from the above linear works, Molteni, Toth 
& Kuznetsov (1999, hereafter MTK) performed 2-D simulations of a shocked adiabatic flow 
and found a non-axisymmetric instability. MTK simulations showed that the instability sat- 
urates at a low level, and a new asymmetric configuration develops, with a deformed shock 
rotating steadily. The mechanism of the instability, however, was not explained in MTK. 
They briefly mentioned a possible link with the non-axisymmetric disc instabilities studied 
by Blaes k Hawley (1988). 

In this paper, we investigate the stability of a shocked accretion flow in the adiabatic 
and inviscid limit. We compare our linear results with the non-linear simulation results of 
MTK and manage to explain the different behaviours in their simulations. The paper is 
organized as follows. In Section 2, we present a set of linearized equations and boundary 
conditions. Subsequently, in section 3, we present the linear results by numerically solving 
the equations. Finally, in section 4, we summarize our conclusions and present a discussion. 
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2 EQUATIONS 

An adiabatic flow around a black hole is considered in the pseudo-Newtonian potential 
introduced by Paczyhski and Wiita (1980), $ = —GM/(r — r g ). Equations are made di- 
mensionless by using the Schwarzschild radius and the speed of light as reference units, i.e., 
r g = 1 and c = 1. In this paper, the thickness of the flow is approximated as a constant for 
the sake of simplicity, as in Nakayama (1994), Blaes (1987), MTK and Paper I. 

The stationary flow is described by the conservation of mass and the Bernoulli equation: 

prv r = const. , (1) 

e = l + ^ + ^T-i(^T) =const " (2) 

where p is the density, v r is the radial velocity, c s = ('jp/p) 1 ^ 2 is the sound speed, / is the 
specific angular momentum, and e is the Bernoulli constant. The structure of a stationary 
flow involving a standing shock can be obtained from the above equations by a given pair 
of (/, e) (see Appendix A of MTK for details). 

The continuity equation and the Euler equation are written as follows: 

^ + v-(H = o, (3) 



dv „ 
— + w x v + V 
at 



v 2 c 2 



^ > ( 4 ) 

7 



7-1 2(r-l)_ 
where w is the vorticity, and S is the entropy. 

Apart from vorticity perturbations and acoustic waves in an isothermal flow, entropy 
perturbations should appear in an adiabatic flow. Thus the linearized equations here are a 
little more complicated than those in Paper I. In order to write out the linearized equations 

in the simplest form, the two functions /, g are defined as follows: 
2 

/ = v r 5v r H —c s 5c s + v^Vp , (5) 

_ 5p 5v r 

9 = — + — , 6 

p v r 

where / is the perturbation of the Bernoulli constant and g is the perturbation of the mass 
accretion rate. The frequency uj' measured in the rotating frame is defined as: 

uj' = uj — mVL , (7) 

where u is the complex frequency of the perturbation, m is the azimuthal wave number, and 
Q = l/r 2 is the angular velocity. With the standard method of linear stability analysis, i.e., 
assuming perturbations to be proportional to e - tuJt + tmi P ^ the following differential system is 
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obtained: 
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where B = rv r w z — imc 2 SS/^f : w z is the vorticity along the rotation axis, M. = —v r /c s is 
the radial Mach number, and the subscript "sh" denotes the shock position. The boundary- 
conditions corresponding to a perturbed shock velocity Av r are obtained: 



/ sh = (V+ - vJ)Av r , 

J_ (\_ V 

W U+ V. 



9sh 
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dr 



(10) 
(11) 

(12) 
(13) 



B s h = , 

where the subscripts "-" and "+" denote the pre-shock and post-shock values, respectively. 
In addition to the two boundary conditions Eqs. (10-11) at the shock, a third equation is 
obtained from the critical condition at the sonic point, 

Uj' f rson i^-dr 

^/so„ + £so„ + SS sh e J ^ «r = , (14) 

where the subscript "son" denotes the sonic point. These three boundary conditions are used 
to numerically solve the differential system Eqs. (8-9) and to determine the eigenfrequencies 

UJ. 

The methods for obtaining the above linearized equations and boundary conditions are 
similar to those in Paper I (see Appendices B and C of Paper I for details). 



3 NUMERICAL RESULTS 

The standard Runge-Kutta method is used to integrate differential equations from the sonic 
point to the shock. In our calculations, the adiabatic index 7 is fixed to be 4/3 as in MTK. 
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Table 1. A sample of 15 shocked accretion flows for numerical calculations. 



case 




I 


M eh 


min 


to = 1 linear results 


MTK simulation results 


\ 


5.3 


1.7770 


0.651 


0.650 


Lllls LaUJt 


VOtTlll^Y* f~\Gf"*lllci"i"1f~\t"*k 


2 


5.3 


1.8000 


0.574 


0.572 


Llllb LcLUJL 


TOfVlll^Y* /~\Gf"*lllf*"i"1/~\t"*k 


3 


5.3 


1.8100 


0.539 


0.536 


unstable 


beating 


4 


5.3 


1.8225 


0.493 


0.489 


unstable 


irregular 


5 


7.8 


1.8000 


0.629 


0.556 


unstable 


regular oscillation 


6 


7.8 


1.8100 


0.595 


0.520 


unstable 


beating 


7 


7.8 


1.8200 


0.560 


0.482 


unstable 


beating 


8 


12.7 


1.8200 


0.684 


0.465 


unstable 


regular oscillation 


9 


17.2 


1.8255 


0.760 


0.439 


unstable 


nearly stable 


10 


23.4 


1.8620 


0.710 


0.294 


unstable 


regular oscillation 


11 


23.4 


1.8720 


0.664 


0.255 


unstable 


leaves domain 


12 


19.9 


1.9200 


0.372 


0.074 


unstable 




13 


34.8 


1.9400 


0.405 


0.034 


unstable 




14 


16.8 


1.8900 


0.481 


0.183 


unstable 




15 


19.8 


1.8000 


0.923 


0.535 


stable 





Thus the lowest post-shock Mach number (corresponding to extremely strong shocks) is: 
•Msh = [(7— 1)/27]5 = 0.354. Tab. 1 presents a sample of 15 shocked accretion flows, of which 
cases 1-11 were exactly taken from Table 1 of MTK. As shown in Tab. 1, for cases 1-11, the 
ranges of the post-shock Mach number A4 s h(0. 493, 0.760) and the minimum Mach number 
in the subsonic region between the sonic point and the shock .M mm (0.255, 0.650) are a little 
narrow compared with the theoretical ones, .M s h(0.354, 1) and .M mm (0, 1), respectively. Thus 
cases 12-15 are added into Tab. 1 to make the ranges wider, i.e. .M s h(0.372, 0.923) and 
A^min(0. 034, 0.650). Our linear results for m = 1 and MTK simulation results are listed on 
the sixth and seventh columns of Tab. 1, respectively. 

Fig. 1 shows the domain (/, e) of angular momentum and Bernoulli constant for which 
an adiabatic flow, subsonic far from the accretor, may be accreted onto a black-hole through 
a stationary shock. Each point between the two solid lines represents an inner shock (post- 
shock acceleration) and an outer shock (post-shock deceleration). Since the inner shock was 
already found to be unstable to axisymmetric perturbations, we will concentrate on the 
stability of the outer shock against non-axisymmetric perturbations. 



3.1 Linear results compared with MTK simulation results 

As shown in Tab. 1, our linear numerical calculations find that cases 1-14 are unstable, 
and case 15 is stable, whereas the non-linear simulations of MTK found that cases 1-11 are 
unstable except for case 9. The agreement of linear and non-linear results can be understood 
as follows. In their simulations, -M sn of the stable flow (case 9) is the largest one (among 
cases 1-11). Similarly, in our linear calculations, A4 s h of the stable flow (case 15) is also 
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I 



Figure 1. The two thick solid lines are the threshold for shock-included solutions. The dotted lines measure the shock strength 
by the value of Ai s h indicated on the left. The thin solid lines correspond to the value of the minimum Mach number Ai m in 
indicated on the right. The 10 crosses represent cases 1-8 and 10-11, which are unstable in MTK simulations. The filled triangle 
represents case 9, which is stable in MTK simulations. The 3 filled squares represent cases 12-14, which together with cases 
1-11 arc linearly unstable. The filled circle represents case 15, which is linearly stable. 

the largest one (among cases 1-15). Thus both linear and non-linear results indicate that 
the shock is generally unstable to non-axisymmetric perturbations except for A4 s h — > 1, i.e. 
very weak shock. We therefore argue that there is no essential difference between linear and 
non- linear results. This general instability is even clear from Fig. 1, which shows that both 
the linearly stable shock (case 15, filled circle) and the non-linearly stable shock (case 9, 
filled triangle) locate very close to the right border, whereas the other unstable shocks stand 
everywhere except for a very narrow region close to the right Ai s h ~ 1 border. 

As a typical example, the eigenspectrum of case 5 is shown in Fig. 2 for perturbations 
< m < 3. The shock is linearly stable to axisymmetric perturbations (m — 0), which coin- 
cides with the criterion "post-shock acceleration causes instability". Despite the post-shock 
deceleration, however, the shock is found to be unstable to non-axisymmetric perturbations 
with m — 1, 2, 3. The MTK simulations found that, for case 5, the perturbed axisymmetric 
shock will finally change into an m = 1 deformed shock. Such a non-linear evolution can 
be well explained by our numerically obtained eigenmodes. As shown in Fig. 2, the fastest 
growth rate correponds to m — 1, u — 0.0725 + 0.00715?, which indicates that the m — 1 
perturbations should dominate over the others. 

To understand the different behaviours of perturbed shocks in MTK simulations, we focus 
on cases 1-7 since there exists continuous change from 1 to 4 and from 5 to 7, respectively. 
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Figure 2. Case 5, eigenspectrum for < m < 3. 
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Figure 3. Cases 1-7, growth rates for 1 < m < 3. 



Fig. 3 shows that the fastest growth rate correponds torn = 1 for cases 1 and 5. For the other 
five cases, however, the fastest growth rate does not correpond to m — 1. In other words, 
among cases 1-7, m = 1 perturbations dominate for cases 1 and 5. As shown in Tab. 1, MTK 
simulations found that the perturbed shock will change into an m = 1 deformed shock for 
cases 1,2 and 5. Comparing the linear results with simulations, we therefore conclude that the 
non-linear behaviours are determined by the fastest growth rate. The m = 1 deformed shock, 
i.e. "regular oscillation", should come into being when the fastest growth rate correponds 
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Figure 4. Case 12, growth rates compared to r ac and T a dv 



to m — 1. On the contrary, other types such as "beating" should appear when the fastest 
growth rate does not correpond to m — 1. 

The only inconsistent solution is case 2. As shown in Tab. 1 and Fig. 3, from case 1 
to 4, the transition direction of the instability is identical for linear and non-linear results, 
i.e., from m — 1 dominance to m / 1 dominance. The transition point locates between 
case 1 and 2 in our linear calculations, however, between case 2 and 3 in MTK simulations. 
As mentioned in MTK, the average distance of the final deformed shock will be a little 
larger than before. Thus the outmoving of the shock may account for the above quantitative 
difference. 



3.2 Instability mechanism 

In the isothermal work of Paper I, the mechanism of the instability was already found to 
be based on the cycle of acoustic waves between the corotation radius and the shock. Such 
a mechanism can be confirmed by the new evidence of Fig. 4, which shows 10 unstable 
eigenmodes for m — 1 of case 12. We choose case 12 because the shock is far away from the 
sonic point and -M mm is low enough, thus the flow has a number of unstable eigenmodes. 
Different from the rough timescales r ac and r ac j v in Paper I, here we numerically calculate 
the exact values of the time of the purely acoustic cycle r ac and the advective-acoustic cycle 
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where the corotation radius r co of the perturbation is defined by uj = mQ , with Q = fi(r co ): 




Fig. 4 shows that the growth time is always around r ac , i.e. c^T ac rs 1, but can be much 
shorter than r a d v , i.e. c^Tadv 1 , which strongly indicates that the mechanism is based on 
the purely acoustic cycle, not the advective-acoustic cycle. However, such an evidence was 
not noticed in Paper I. 

To have a global view, for m — 1, the most unstable modes of all the 15 flows are included 
in Fig. 5, which shows that for low .M m i n , i.e. r ac <C r adv , the most unstable mode well 
matches c^r ac ps 1. This result again confirms the above mentioned mechanism in adiabatic 
flows. The value of WjT ac , however, evidently decreases as .Mmin — ► 1, which implies that the 
instability is suppressed by advection. Different from Fig. 5 of Paper I, in Fig. 5 here we 
choose .M m in as abscissa instead of -M sn because an adiabatic flow with 7 = 4/3 can not 
have very low -M sn , and more importantly, -M min implies the strength of advection. 
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Figure 6. The positions of the corotation radius of the most unstable mode are indicated for 15 adiabatic flows (circles) and 
24 isothermal flows (squares). The dotted line corresponds to (r co — »*son)/(r s h — r BO n) = 1 — -M s h- 

3.3 Oscillation period 

The growth rate is relevant to the imaginary part of the eigenfrequency u>i, whereas the 
oscillation period is relevant to the real part uo r (or r co ). Fig. 6 includes all the 15 flows in 
Tab. 1 together with the 24 isothermal flows from Paper I, which shows the relationship 
among r co , r son , r sh and M.^. For M.^ that is not very low, Fig. 6 suggests the following 
good approximation: 

^ co '"son -\ \ a / -i o \ 

= 1-M sh . (18) 

7"sh ''"son 

Thus an analytical formula for calculating the oscillation period is obtained: 



p _ 2 ^co 



2tt 



Fsh 



(?"sh - r son)-M sh ] 



(19) 



I I 

The lowest post-shock Mach number A4 s h = [(7 — l)/2^ implies that only an isothermal 
flow or a flow with 7 — > 1 can have extremely low value of M s h- Thus Eq. (19) should work 
well in normal adiabatic flows with 4/3 < 7 < 5/3. In addition, we define the linear period 
Pun for m — 1 as follows, 
2tt 



lin 



(20) 



where u;™ ax is the real part of the eigenfrequency corresponding to the most unstable mode. 

The different types of periods are shown in Fig. 7. The lower and the upper dotted 
lines correspond to the rotation period at the sonic point P son = 2irr1 on /l and at the shock 
P sn = 2nr^ h /l, respectively. This figure shows that both the linear Pi in and the non-linear 
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Figure 7. Three types of periods for cases 1-15: the analytical .Pan; the linear Pn n and the MTK simulation results -Psim- 

The 

upper and the lower dotted lines correspond to P s h and P S on, respectively. 



P S im always locate between P son and P sn (except for case 2), which is in good agreement with 
the condition that the corotation radius should locate between the sonic point and the shock. 
Fig. 7 demonstrates that the analytical period P an is a good approximation for the linear 
period Pu n and not far from the non-linear period P snn . As mentioned in MTK, the average 
distance of the final deformed shock will be a little larger than before. For example, in case 
5, the final distance r sh varies between 9-11, which is larger than the original value of 7.8. 
If this outmoving is taken into consideration, P sn = 2nr^ h /l should therefore be even larger 
thus P sim of case 2 should indeed locate between P son and P sn . Futhermore, the outmoving 
well explains why P sim is always larger than p in and P an (as shown in Fig. 7). 



4 CONCLUSIONS AND DISCUSSION 

The main results can be summarized as follows: 

(i) The general non-axisymmetric instability of the outer shock, which was previously 
found by non-linear simulations for adiabatic flows and linear calculations for isothermal 
flows, is confirmed in the present paper by linear calculations for adiabatic flows. 

(ii) The simulation results of MTK are well explained by our numerically obtained linear 
eigenmodes. 

(iii) New evidence is shown to support the argument in Paper I for the mechanism that 
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the instability is based on the cycle of acoustic waves between the corotation radius and the 
shock. 

(iv) An analytical formula is obtained to calculate the oscillation period P from the 
physical parameters r son , r sh , / and A4 s h- 

The present work is for an inviscid accretion flow with a standing radial shock. Realistic 
astrophysical accretion flows must have viscosity, and it remains controversial whether a 
standing shock can indeed form in a viscous accretion flow around a black hole. In an 
accretion flow, the gravitational potential energy is converted to kinetic and thermal energy 
of the accreting gas. Based on the energy consideration, Narayan, Mahadevan & Quataert 
(1998) identified three regimes of accretion: the radiative cooling-dominated accretion flow, 
the advection-dominated accretion flow (ADAF), and the low energy generated accretion 
flow. For cooling-dominated flows or ADAFs, the released gravitational potential energy is 
mainly converted to thermal energy by viscous stresses and then radiated away or advected 
into the central black hole. Except for the region very close to the black hole (<5r g ), the 
radial motion of these two kinds of accretion flows keeps subsonic because of the low kinetic 
energy. This is the physical reason why there are no shocks in the ADAF-thin disc solutions 
(e.g. Narayan, Kato & Honma 1997, hereafter NKH), and this result does not rely on the 
mathematical technology. In fact, Lu, Gu & Yuan (1999) recovered the shock-free ADAF-thin 
disc solutions using the Runge-Kutta method (by adjusting the specific angular momentum 
j accreted by the black hole and the sonic point r son of an ADAF to match a thin disc at 
the outer boundary), which are identical with those in NKH using the relaxation method 
(by assuming the value of the outer boundary r out and the conditions at r out , and then 
solving the set of euqations and calculating out j and r son as eigenvalues). In our opinion, a 
standing shock may form if the original flow belongs to the third regime mentioned above, 
i.e. a low energy generated accretion flow. If the flow has very low angular momentum / 
at large distance, i.e. close to Bondi accretion (in an ADAF I is not very low, see Fig.2 of 
NKH), the gravitational force dominates over the centrifugal force. The flow is accelerated 
efficiently in the radial direction and becomes supersonic far from the central black hole 
(Yuan 1999), then a shock is likely to develop. Such a flow is a low energy generated one 
because the released gravitational potential energy is mainly converted to kinetic energy 
rather than thermal energy. NKH has also agreed that this is the situation in which a shock 
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can indeed be physical. Thus we believe that flows with low angular momentum and weak 
viscosity (i.e. low energy generation) will have shocks. 

It is well-known that the PPI occurs with either an inner or an outer reflecting boundary, 
or even more efficiently with both. The boundary in the present paper is a standing shock. 
However, a shock is not the only type of boundary which can reflect acoustic waves. For 
example, the transition from a Shakura-Sunyaev disc (Shakura & Sunyaev 1973) to an ADAF 
ought to be very sharp (NKH, Manmoto & Kato 2000, Lu, Lin & Gu 2004). Such a sharp 
transition surface can also reflect acoustic waves, thus the PPI may occur and result in the 
QPO in this system. We therefore argue that the QPO should be a common phenomenon 
in accretion flows with angular momentum. In particular, the QPO frequencies in a 3:2 
ratio in black hole X-ray binaries (McClintock & Remillard 2005) may be explained by this 
instability for both m=2 and m=3 dominance, such as cases 2 and 7 (as shown in Fig. 3). 
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